Empirical Bayes single nucleotide variant-calling for next-generation sequencing data

One of the fundamental computational problems in cancer genomics is the identification of single nucleotide variants (SNVs) from DNA sequencing data. Many statistical models and software implementations for SNV calling have been developed in the literature, yet, they still disagree widely on real datasets. Based on an empirical Bayesian approach, we introduce a local false discovery rate (LFDR) estimator for germline SNV calling. Our approach learns model parameters without prior information, and simultaneously accounts for information across all sites in the genomic regions of interest. We also propose another LFDR-based algorithm that reliably prioritizes a given list of mutations called by any other variant-calling algorithm. We use a suite of gold-standard cell line data to compare our LFDR approach against a collection of widely used, state of the art programs. We find that our LFDR approach approximately matches or exceeds the performance of all of these programs, despite some very large differences among them. Furthermore, when prioritizing other algorithms’ calls by our LFDR score, we find that by manipulating the type I-type II tradeoff we can select subsets of variant calls with minimal loss of sensitivity but dramatic increases in precision.

data, the type of mutations sought (somatic or germline), noise characteristics of the dataset, etc.As such, most users/analysts may not have enough expertise to adjust those pre-determined parameters properly.On the other hand, relying on default values for parameters may result in unreliable results and is problematic when comparing programs, because their default parameters differ.
The above approaches are either conventional Bayesian or frequentist methods, and do not take multiplicity and testing efficiency issues into account.Moreover, many programs fail to output a criterion by which SNV calls can be ranked, such as p-values, so that adjusting preference between type I and type II errors or constructing ROC curves is not naturally supported.
Some studies in the literature proposed germline variant detection in multi-sample NGS data using Bayesian methods.For example, Huang et al. 11 raised the question of how to quantify sequencing errors.They argue that it may not be ideal to quantify sequencing error only based on quality scores and propose an empirical Bayes method based on a Binomial model to estimate sequencing errors and identify variants.In another study, Zhao et.al 12 developed an optimal empirical Bayesian testing procedure to detect variants in NGS data, which is based on pooling a normalized amount of DNA from multiple samples.Due to a capacity issue, they assume samples may be distributed and sequenced independently in M > 1 pools and each pool consists of N individu- als.Although having M pools with N individuals may be cost-effective, however, our available samples (similar to many other clinical labs) instruct us to base the model only on one DNA sample for one individual.Similar to Huang et al. 11 , Ding et al. 13 suggested that estimating sequencing errors based on quality scores may not be precise.Using single-sample NGS data, they proposed an empirical Bayes method to model the dispersion of minor AF and call germline variants.
In this paper, we propose a Binomial model and include an error rate based on a weighted average of base call and mapping quality scores as part of model parameters.We use an empirical Bayesian approach to develop a local false discovery rate (LFDR) estimator for SNV detection.In contrast to existing algorithms, our novel approach calls SNVs not just on a site-by-site basis, but by simultaneously using information across all the sites to build a probabilistic model of the data.
To the best of our knowledge, LFDR estimation has not been employed in the task of variant-calling, but it has been well developed in a variety of other contexts, and different strategies have been introduced in the literature.To name a few LFDR estimation based approaches, readers may refer to Pan et al. 14 , Efron et.al 15 , Efron 16 , Padilla and Bickel 17 , Yang et al. 18 , Karimnezhad and Bickel 19 , and Karimnezhad 20 .
Obviously, all variant-calling algorithms determine mutations based on some observed evidence, but since evidence supporting variants (including number of reference and alternative read counts) differs from site to site, the confidence in calling a site as a variant site varies.We fill this gap by introducing an LFDR-based algorithm that meaningfully scores variants called by any variant caller and prioritizes them from most to least probable variants.This helps with significantly reducing many false positives.
The structure of this work is as follows.In the "Methods" section, we provide a detailed presentation of the model, method and the algorithm we propose for germline SNV calling.We also propose a modified version of the LFDR algorithm so that it can prioritize variants called by any desired variant caller.In the "Results" section, we use a suite of gold-standard cell line data to evaluate the performance of our proposed LFDR approach against a collection of widely used, state of the art programs: MuTect2, SAMtools, VarScan2, VarDict, and Pisces.

Methods
We consider analyzing a single DNA sample from a healthy patient, where only germline mutations are expected to be seen.We also assume that the data has been mapped to a reference genome, which specifies one of A, C, G, or T bases for every position.Each mapped read shows one of the four possible bases (A, C, G, or T) for the same position.With no loss of generality, we also suppose, for a replicate (technical or biological), that sequencing covers p sites (individual positions in the genome) of which p 0 sites are non-mutant.For each locus i, i = 1, . . ., p , we assume that there are K i (short) known number of reads covering the locus i.At each locus i, we suppose there are four possible bases (A, C, G, T) of which R i random reads carry the reference allele.We further suppose that of the remaining K i − R i reads, M i random number of reads carry the most alternative (dominant) allele and, we assign the remaining random alternative reads to X 1i and X 2i .We briefly refer to these random variables by X i = (R i , M i , X 1i , X 2i ) .Generally, the larger M i is, the more support we have for an alternative allele being present at site i in the DNA being sequenced.Indeed, if the data had no sequencing or mapping errors, M i > 0 could only arise as the result of an alternative allele being present.However, it is well-recognized that sequence data does contain errors.Although those errors can arise in many ways, the standard model error employed by the community is the error in mistakenly reporting an alternative allele as one of the other three possible alleles, and vice versa.We refer to this error by e i .Obviously, if there is no error, mutant sites are expected to be those p − p 0 sites, for which M i is positive.Suppose that θ i represents the dominant alternative AF at the ith site.Obviously, given the information that ith site has no mutation, θ i = 0 and otherwise, θ i > 0 .We assume P(θ i = 0) = π 0 and P(θ i > 0) = 1 − π 0 .If θ i > 0 , we assume it follows a probability density function g(•) .We assume that each read is drawn independently from either the reference or alternative allele pools.However, we also assume an independent chance e i of read error in each read at each site, which randomly changes the correct allele to one of the other three options.The error term e i can be calculated in different ways.For example, in Mutect, error is estimated based on Phred- scaled quality scores of mapped reads to each position.We suggest to take the error rate as a weighted average of average base call qualities and average mapping qualities.Supposed at position i, for an arbitrary base pair U ∈ {A, C, G, T} , BQ U i and MQ U i represent average base call quality and average mapping quality across all reads at that position, respectively, and let Q U i be the average of BQ U i and MQ U i , i.e., . Depend- ing on whether the base pair observed at position i is a reference (R), dominant (M) or one of the remaining alternative alleles ( i , and use them to calculate the weighted average quality . Then, we define the error rate to be e i = 10 − Q i 10 .As in Mutect 6 , we define the following four probabilities.The probability that a randomly chosen read covering site i shows the reference allele is . This is explained as the sum of the probabilities that the DNA actually contained the reference allele at site i and it was correctly read, (1 − θ i )(1 − e i ) , and probability that the DNA actually contained the alternative allele but it was misread, and by chance, it was misread as the reference allele, θ i e i 3 .Similarly, the probability that a random read covering site i shows the alternative allele is 3 , which arises either as correct reading of the alternative allele or misreading the reference allele as the alternative.Finally p X 1i = p X 2i = e i 3 is the chance that one of the two other alleles that are not reference or alternative occurs at site i.Putting these together, the probability of the total data at site i, X i is multinomial with K i tries (reads) and probabilities of reference, alternative, and the other two alleles.At this point, we diverge from the model of Mutect by proposing a new model governing the parameter θ i , which ultimately allows us to develop our empirical Bayesian LFDR approach.We assume , where θ i could be either 0 or positive with probabilities π 0 and 1 − π 0 , respectively.If θ i > 0 , we assume θ i ∼ g(•) , where g(•) is unknown.In our proposed approach we do not need to know the form of g(•).
To discover whether site i is a mutant site, we propose testing the null hypothesis H 0i : θ i = 0 against the alternative hypothesis H 1i : θ i > 0 , i = 1, . . ., p .To do so, we focus on estimating ψ i ≡ P(θ i |X i ) , the posterior probability that the null hypothesis H 0i is true.Once estimated, it is compared with a pre-specified threshold leading to either rejecting or failing to reject the null hypothesis.The quantity ψ i is well-known in the literature as LFDR.
With the above settings, the following probability functions are derived under the null and the alternative hypotheses, respectively, and Now, using the Bayes formula, the LFDR can be expressed by ( 1) (3) where P(X i |θ i = 0) and P(X i |θ i > 0) are given by ( 1) and ( 2), respectively, and π 0 is the proportion of non- mutant sites.Both the parameters π 0 and g(θ i ) (in P(X i |θ i > 0) ) are unknown and need to be estimated before making any inference.Now, let δ i be a binary decision rule corresponding to ith set of hypotheses H 0i and H 1i .We assume that δ i = 1 if the null hypothesis H 0i is rejected, and δ i = 0 , otherwise.But since such binary decisions can lead to some errors, we define the following loss function when testing the null hypothesis at site i: where l I and l II are loss values incurred due to making type I and type II errors, respectively, see Karimnezhad and Bickel 19 .Also, see Zhao et.al 12 where a specific version of this loss with l I = and l II = 1 is used.
To take all p sites in a sample to account, suppose that δ = (δ 1 , . . ., δ p ) represents a vector of estimators of θ = (θ 1 , . . ., θ p ) and measure the aggregated loss by L(θ , δ) = p i=1 L(θ i , δ i ) , which takes both type-I and type-II errors into account.Now, to derive a Bayesian decision rule in the above-mentioned hypothesis testing problem, let ρ X (θ, δ) = E[L(θ, δ)|X] denote the posterior risk of choosing a decision vector δ , where X = (X 1 , . . ., X p ) .Then, similar to Karimnezhad and Bickel 19 , it can be verified that the Bayesian decision vector where ψ i is given by (3), minimizes the posterior risk w.r.t.δ.
The above equation suggests that mutant sites can be determined by estimating ψ i , i = 1, . . ., p , and then comparing it with the threshold l I l I +l II .However, estimating ψ i is challenging due to the complicatedness of the form of P(X i |θ i > 0) as well as g(•) , the distribution of AFs.
One immediate solution to estimating g(•) would be to assign a non-informative prior to the alternative AFs θ i .In a similar situation, Zhao et.al 12 consider a uniform distribution on the interval (0, a), and then they estimate the value of a as well as π 0 using the traditional method of moments.While simple, this assumption seems unrealistic.Figure 1 represents the distribution of AFs at known mutant as well as non-mutant sites from TST-GM12878 data (see Results for more information).
From Fig. 1, we observe that the assumption that AFs take any number between 0 and 1 with equal probability is unrealistic.Estimating π 0 may be connected to the estimation of g(•) and consequently, a bad estimate of g(•) may directly lead to a bad estimate of π 0 .This indeed is true in some datasets with ambiguous AF ranges, including our datasets.Clearly, from this figure, AFs of mutant sites is distributed around either 0.5 or 1.This highlights the importance of finding a suitable estimate of g(•) .Dependency of estimated values of parameters of LFDR to each other can be seen in different algorithms in the literature.See for example Zhao et.al 12 where their proposed estimator of π 0 is a function of the estimated value of the upper bound a in the uniform prior assigned to the alternative AF.For other examples, see Pan et al. 14 and Karimnezhad and Bickel 19 .
As Fig. 1 suggests, one may fit a beta distribution to model the alternative AFs, which adds extra parameters to the model.To avoid model complexity, as an innovative approach, we suggest estimating g(•) empirically, An immediate solution to the above problem would be to take all sites for which M i counts are positive.But as many non-zero alternative AFs in NGS studies turn out to correspond to artifacts, we need to find and exclude them from our list so that a precise empirical distribution for g(•) can be estimated.Suppose that I s stands for the set of s indices of those non-zero θ i corresponding to mutant sites, and let θ s be the vector of the correspond- ing empirical AFs.In our settings, s represents the number of elements of the set I s .Now, let G(•) represent the cumulative distribution function (CDF) that corresponds to g(•) .Then, the corresponding empirical CDF is given by G s (t) = 1 s l∈I s 1 θ l ≤t , t ∈ ℜ .Obviously, due to the strong law of large numbers, G s (t) converges to G(t) as n → ∞ almost surely, for every value of t, and this implies that G s (t) is a consistent estimator.Then, one may estimate P(X i |θ i > 0) by Since I s includes many artifacts, s needs to be optimally learned.We propose an algorithm that allows for repeat- edly updating the set of mutant sites and leads to a precise estimate of g(•).

An empirical Bayes mutation detection procedure
To discover mutant sites, we propose the procedure outlined in Algorithm 1.The proposed algorithm, like many existing procedures in the literature, including Mutect2, VarScan2 4 , VarDict 7 and Pisces 8 is comprised of three main steps.
In the pre-processing step, we propose excluding low quality bases from the input file (either FASTQ or BAM format).Thus, observations in our model are only those bases that pass a minimum average BQ and average MQ threshold.BQ corresponds to an error rate of e = 10 − BQ 10 , see for example Cibulskis et al. 6 and Dunn, et al. 8 .Thus taking BQ = 20 corresponds to an error rate of 0.01 which can be inferred as expecting one miscalled base in reading 100 bases.MQ is related to aligning reads to a reference genome, and it usually varies between 0 and 60.Similar to Pisces, we impose a threshold of 20, and 30 to the minimum BQ and MQ, respectively.Different publicly available softwares such as SAMtools 3 and bam-readcount (https:// github.com/ genome/ bam-readc ount) can be applied to exclude such low quality reads.Once low quality bases are excluded, the counts K i , R i , M i , X 1i and X 2i need to be formed.
Although in Eq. ( 5) we only include all non-zero empirical AFs, one may tweak the algorithm by just assuming that AFs follow, for example, a uniform distribution on the interval (0, 1), or a subset of it.Thus, I j s in Step 2.6 of Algorithm 1, may be alternatively taken to be a set of N (for example 1000) sites for which their AFs are randomly generated from a uniform distribution on the interval (0, 1).However, once the estimation procedure enters Step 2.10, the updated I j s set will only depend on the actual empirical AFs.A third approach would be to estimate LFDRs by focusing only on the original uniformly sampled AFs, without updating the set of indices I j s in step 2.10, i.e., for all j, I j s = I 1 s .We refer to this approach by "uniform" estimation.For comparison purposes, we use these three methods in our data analysis and will refer to them by "empirical", "uniform/empirical" and "uniform" estimation of g(•) , respectively.In Step 2.10, we took ǫ = 0.001.
In the post-processing step, we take 0.01 and 10 as default AF threshold (AFT) and read depth threshold (DPT) values, respectively.This is not uncommon, and many variant callers impose such thresholds.For example, Pisces 8 has the same default values for these thresholds.We remark that one may apply this filtering in Step 1.However, the final estimated π 0 will reflect the proportion of non-mutant sites w.r.t. the filtered input file rather than the original one. (5) X1i+X2i .

LFDR as a variant prioritization tool
We suggest that our LFDR calculation can be applied to prioritize any list of mutations called by any mutation caller from most to least probable mutations.This can be done by modifying Algorithm 1 so that π 0 is estimated as the number of reference sites according to the caller, out of all those assayed, and g(•) is estimated as the empiri- cal distribution of variant sites.The LFDR calculation is then applied to all variant sites using these estimates of π 0 and g(•) .Algorithm 2 outlines the prioritization steps. Vol

Results
We now focus on evaluating our proposed mutation calling algorithm on some data generated by clinical assays.We employ DNA sequencing data generated from two Coriell cell lines, GM12877 and GM12878, previously studied by Karimnezhad et al. 9 .These are immortalized healthy (non-cancer) cell lines that are expected to carry only germline variants.These cell lines have been extremely well characterized by numerous groups.Thus, we can know exactly which variants they harbor relative to the reference human genome, allowing us to accurately assess true and false positives and negatives.In particular, we rely on a list of known mutations published by Eberle et al. 22 where a comprehensive and genome-wide catalog of high-confidence variants mutations for a collection of Coriell cell lines, including GM12877 and GM12878 is presented (available through https:// www.ncbi.nlm.nih.gov/ gap/ under accession phs001224.v1.p1).The data we used was sequenced in two different ways: (1) on an Illumina NextSeq500 sequencer using Illumina's TruSight170 targeted sample preparation method (TST170 for short), and on an Ion Torrent PGM sequencer using an Oncomine Focus targeted panel (OF for short).The TST170 data spans 514761 total bases, covering parts of 170 genes, and includes six technical replicates of each of the Coriell cell lines GM12877 and GM12878.By intersecting the list of known mutations with the TST170 genomic regions, we determined that our GM12877 and GM12878 data should contain 336 and 343 mutations, respectively.The OF panel covers 29008 total bases in 47 genes, and includes three technical replicates of each of the Coriell cell lines GM12877 and GM12878.There should be 24 and 26 known mutations present in each replicate of the GM12877 and GM12878 data, respectively.For more information regarding the datasets and also sequencing platforms, readers may refer to Karimnezhad et al. 9 .
For a given single replicate, denote the total number of true positives, false positives and false negatives by TP, FP and FN respectively.We measure the performance of the algorithms by computing precision or positive predictive value Prec = TP TP + FP and sensitivity Sens = TP TP + FN .A good algorithm is expected to have high Prec and Sens rates.
To apply the algorithm on the TST170 as well as OF replicates, we used the bam-readcount package to calculate the counts K i , R i , M i , X 1i and X 2i , and then followed Steps 1-3 of Algorithm 1.We then compared the list of final detected variants with the lists of known mutations.Also, to investigate the impact of a chosen LFDR threshold in Step 2.1 of the algorithm on the detection accuracy, we picked 10 −300 , 10 −200 , 10 −100 , 10 −50 , 10 −20 , 10 −10 , 10 −2 and 0.5.
We measured the performance of the proposed approaches by calculating TP, FP, FN, and consequently Prec and Sens values for different datasets and LFDR thresholds.As a second performance measurement, we compared the proposed algorithm with Mutect2, VarScan2, SAMtools, VarDict and Pisces.Illumina offers Pisces for the analysis of TST170 data through their customized pipeline.The application accepts paired-end fastq files as inputs, generates BAM files, and after aligning to the reference human genome (hg19) by the Isaac aligner 23 , uses Pisces to generate a list of mutations.Then, the final list of mutations is generated after some internal filtering.Our work was carried out in the context of a clinical sequencing project, where use of hg19 remains widespread, prompting us to use hg19 for alignments.However, we would expect very similar results with hg38.Our data is obtained by targeted sequencing of a panel of clinically-relevant genes, and there should be little or no difference between hg19 and hg38 for these highly-studied genes.We should add that, to compare the performance of our proposed algorithm with Pisces, as well as VarScan2, Mutect2, SAMtools and VarDict, we used the same aligned BAM files to reduce some possible alignment errors.Because differences in default parameters of the five algorithms in Table 1 are potential sources of discrepancies in the list of final mutations, we set their parameters as similar as possible.We set the minimum variant AF to 0.01.We set minimum base call quality and minimum mapping quality to 20.Finally, we set the minimum coverage for a called variant at 10 reads.
Figure 2 represents the average Prec and Sens values over replicates based on calls made by the LFDR approach for different thresholds, as well as the other five variant callers.From this figure, we observe that for all the LFDR thresholds, Sens values are mostly high while Prec is increased by decreasing the LFDR threshold.Indeed, when the LFDR threshold is high, say 0.5 for example, many false positives are allowed to be in the list of variants called.However, a strict LFDR threshold, say 10 −50 for example, leads to few(er) false positives and consequently an improved Prec is gained.All the three estimation methods (empirical, uniform/empirical and uniform) performed nearly the same.Especially, when the LFDR threshold is below 10 −50 , the three methods led to almost the same Prec and Sens values.
Comparing with the other variant callers, we observe that the LFDR approach has equal or better Sens for some LFDR thresholds compared to all other algorithms.And for other thresholds, it has equal or better Prec.For example, from panel (a) in Fig. 2 we observe that the LFDR approach with a threshold of 10 −100 outperforms all the other five variant callers in terms of Prec.However, it reports a smaller Sens.When looking at other panels, slightly different but remarkable performance is observed.For example, by choosing the same LFDR threshold, we observe in panel (b) that the LFDR approach outperforms MuTect2 in terms of Prec, and has almost the same Sens, and outperforms the remaining four variant callers in terms of Prec.However, its Sens is a bit smaller than the others.But when looking at panel (d), we observe that the LFDR approach with the same LFDR threshold outperforms all the other variant callers in terms of both Prec and Sens.We observe that having a varying LFDR threshold allows for a wide range of tradeoff between Sens and Prec.
Figure 3 represents convergence of estimated π 0 in replicate 1 of each dataset for different LFDR thresholds.From this figure, we observe that the algorithm converged mostly in 3-7 iterations, and most remarkably, for stringent LFDR thresholds ( 10 −100 or less), it converged in at most 4 iterations.We also note that in each panel, those stringent LFDR thresholds resulted in the greatest accuracy in estimating π 0 .In fact, the lower is the LFDR threshold, the more accurate is the corresponding estimated π 0 .Comparing with Fig. 2, we also observe that the most accurate estimated π 0 leads to the highest Prec.However, this is not neccessarily true when looking at Sens values.We also investigated the impact of each method on the magnitude of LFDR estimates.Figure 4 represents estimated LFDRs for replicate 1 of each of our datasets based on taking 10 −300 as the LFDR threshold.This figure reflects that both the empirical and uniform/empirical approaches coincide on detecting the same variants.The corresponding estimated LFDRs are mostly either close to zero (that corresponds to mutant sites) or one (that corresponds to non-mutant sites), revealing that both approaches were able to nicely categorize the sites as either mutant or non-mutant sites.But, LFDRs estimated using the uniform approach did not follow this structure and there are many sites with medium estimated LFDRs.This concludes that both the empirical and uniform/empirical approaches outperform the uniform approach in classifying sites to non-mutant as well mutant categories.
Following the steps in Algorithm 2, we calculated the LFDR values for variants detected by the five variant callers (MuTect2, VarScan2, SAMtools, VarDict and Pisces) across replicate 1 of all the datasets.Figure 5 represents the corresponding LFDRs and whether the variants are TPs or FPs, along with the number of TPs and FPs.From the figure we observe that all those TPs led to either zero or very close to zero estimated LFDR values, as expected.For example, we observe from panel (a) that of the total 1532 variants detected by MuTect2,  We used three methods to estimate g(•) in Eq. ( 2).As expected, both the empirical and uniform/empirical approaches performed the same and identified the same variants.This is because both approaches lead to the same I s set.The uniform approach also led to very similar results.This could be because we are ignoring the magnitude of LFDR values and we just compare them with the chosen thresholds.One possible explanation could also be the fact that the data we used in our project do not have any genuine mutations at low AFs, and thus, accurately estimating g(•) has little effect on the performance.However, as we showed in Fig. 4, the uniform approach is not able to firmly categorize sites to either mutant and non-mutants sites.
In this paper, focused on germeline SNV calling.However, as the proposed model learns the form of g(•) in Eq. ( 2), the proposed SNV calling approach may be used in a tumor-only mode, without comparison to healthy tissue from the same patient, which is often the case in clinical sequencing.The accuracy of the proposed method for tumore-only samples remains open to be investigated.
One may incorporate uncertainty about the parameter θ i using other Bayesian framewroks.For example, one may compute the Bayes factor BF = P(X i |θ i =0) P(X i |θ i >0) and re-express the LFDR quantity defined in Eq. ( 3) as or  where PO = π 0 1−π 0 can be interpreted as the prior odds of having no mutations, see Wakefield 21 for more details.However, having the BF alone does not seem to be enough to identify variant sites and π 0 still needs to be some- how estimated.The beauty of our proposed approach is that it estimates both P(X i |θ i > 0) (and hence the Bayes factor) and π 0 without having any prior information regarding the distribution of θ i .
We end our discussion with highlighting that in some situations, one may notice conflicts in calling a site a variant site.This may happen when multiple replicates are studied simultaneously and a variant caller detects a mutation at a specific site in one replicate and does not call that variant in the other replicate(s).One simple way to resolve this inconsistency would be to take the intersection of the variants called in all replicates.However, this may exclude some important variants (if borderline) and in general ignores the confidence any one replicate gives us.In an extreme case, if a replicate failed to have data in a region, and thus no variants were called, that lack of data would effectively overrule the positive data in other replicates.Alternatively, one may seek to take advantage of advanced decision-theoretic approaches.In a different but still applicable context, Karimnezhad and Bickel 19 developed an empirical Bayesian approach that takes advantage of prior information from multiple reference classes and leads to a unique decision in conflicting situations.This could be a potential future research problem.

Figure 3 .
Figure 3.Estimated π 0 in replicate 1 of each dataset for different values of LFDR thresholds.The solid grey line represents the true π 0 .

Figure 6 .
Figure 6.Average of Prec and Sens values for calls prioritized by the LFDR approach.The pink symbols refer to the Prec and Sens values of each individual variant caller calculated in Fig. 2.

Table 1 .
Default values for options in different mutation callers.A dash means that the corresponding parameter was not defined in the caller's settings.
based on the alternative AFs at a set of known or suspected mutated sites.This raises the question of which set of sites to use.